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Abstract.  In  reliability  engineering  focused  on  the  design  and  optimization  of  structures,  the 
typical  measure  of  reliability  is  the  probability  of  failure  of  the  structure  or  its  individual  components 
relative  to  specific  limit  states.  However,  the  failure  probability  has  troublesome  properties  that 
raise  several  theoretical,  practical,  and  computational  issues.  This  paper  explains  the  seriousness 
of  these  issues  in  the  context  of  design  optimization  and  goes  on  to  propose  a  new  alternative 
measure,  the  buffered  failure  probability,  which  offers  significant  advantages.  The  buffered  failure 
probability  is  handled  with  relative  ease  in  design  optimization  problems,  accounts  for  the  degree 
of  violation  of  a  performance  threshold,  and  is  more  conservative  than  the  failure  probability. 
Keywords:  Failure  probability;  Structural  Reliability;  Reliability-based  design  optimization. 

1  Introduction 

Civil,  mechanical,  naval,  and  aeronautical  structures  such  as  bridges,  building,  offshore  platforms, 
vehicle  frames,  ship  hulls,  and  aircraft  wings  are  subject  to  uncertain  loads,  environmental  condi¬ 
tions,  material  properties,  and  geometry.  It  is  widely  recognized  that  these  uncertainties  must  be 
accounted  for  in  the  design,  maintenance,  and  retrofit  of  such  structures.  The  theory  of  structural 
reliability,  see,  e.g.,  [5],  provides  an  analytic  framework  for  assessing  the  reliability  of  a  structure  as 
measured  by  its  failure  probability  to  be  defined  precisely  below.  The  failure  probability  is  widely 

promoted  to  designers  and  building  code  developers  as  a  tool  for  assessing  and  comparing  designs 
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and  has  successfully  been  applied  to  many  applications,  see,  e.g.,  [5].  While  the  failure  probability 
is  of  significant  importance,  it  also  possesses  troublesome  properties  that  raise  several  theoretical, 
practical,  and  computational  issues.  In  particular,  these  issues  surface  when  the  failure  probability 
is  used  in  design  optimization  of  structures  and  may  lead  to  poor  numerical  performance  of  stan¬ 
dard  nonlinear  optimization  algorithms  such  as  SNOPT  [9],  LANCELOT  [4],  and  NLPQL  [28].  In 
this  paper,  we  discuss  these  issues  and  propose  an  alternative  measure  of  reliability  that  we  call  the 
buffered  failure  probability.  The  buffered  failure  probability  is  handled  with  relative  ease  in  design 
optimization  problems,  accounts  for  the  degree  of  violation  of  a  performance  threshold,  and  is  more 
conservative  than  the  failure  probability. 

The  failure  probability  and  the  buffered  failure  probability  are  defined  in  terms  a  limit-state 
function  g(x,v)  that  is  a  function  of  a  vector  x  =  (xi,  X2, ...,  Xn)'  of  design  variables  (with  prime 
'  denoting  the  transpose  of  a  vector),  which  may  represent  member  sizes,  material  type  and  qual¬ 
ity,  amount  of  steel  reinforcement,  and  geometric  layout  selected  by  the  designer,  and  a  vector 
V  =  (ui,  ^2, ...,  of  quantities,  which  may  describe  loads,  environmental  conditions,  material 
properties,  and  other  factors  the  designer  cannot  directly  control.  The  quantities  v  are  usually 
subject  to  uncertainty  and  their  values  are  therefore  not  known  a  priori.  The  limit-state  function 
represents  the  performance  of  the  structure  with  respect  to  a  specific  criterion  referred  to  as  a  limit 
state.  As  commonly  done,  we  describe  these  quantities  by  random  variables  V  =  (Li,  V2,  •••,  Pm)^ 
with  a  joint  probability  distribution  which  is  regarded  as  known,  although  it  might  need  to  be  esti¬ 
mated  empirically.  To  distinguish  between  the  random  variables  and  their  realizations,  we  denote 
the  former  by  capital  letters  and  the  latter  by  lower  case  letters.  For  a  given  design  x,  gf(x,  V)  is 
a  random  variable  describing  the  (random)  performance  of  the  structure.  We  refer  to  this  random 
variable  as  the  state  of  the  structure. 

By  convention,  gf(x,  v)  >  0  represents  unsatisfactory  performance  of  the  structure  with  respect 
to  the  limit-state  function  and,  consequently,  the  event  {^(x,  V)  >  0}  is  the  set  of  realizations  of 
the  random  vector  V  corresponding  to  “failure.”  We  refer  to  this  set  as  the  failure  domain.  We 
note  that  failure  may  not  necessarily  imply  total  collapse  of  the  structure,  but  may  simply  mean 
the  violation  of  a  prespecified  threshold  for  crack  width,  deflection,  vibration,  etc. 
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The  current  approach  to  structural  reliability  defines  the  failure  probability  of  a  structure  with 
limit-state  function  g(x,  v)  as  the  probability  that  the  state  of  the  structure  takes  on  a  positive 
value.  As  the  failure  probability  depends  on  the  design  x,  we  denote  it  by  p(x).  That  is, 

p(x)  =  P[g{x,Y)>0] 

=  y  ...y v)  >  0)/v(v)dui...(ium,  (1) 

where  /v(v)  is  the  joint  probability  density  function  for  V  and  I{g{x,v)  >  0)  is  the  indicator 
function  defined  to  be  one  if  g{x,  v)  >  0  and  zero  otherwise. 

Our  definitions  of  unsatisfactory  performance  and  the  failure  domain  deviate  in  two  minor  ways 
from  those  of  some  other  authors,  see,  e.g.,  [5].  First,  we  exclude  the  realizations  v  corresponding  to 
g{x,  v)  =  0  from  the  failure  domain.  Of  course,  if  the  probability  of  the  event  {(/(x,  V)  =  0}  is  zero, 
as  is  typically  the  case  when  V  are  continuous  random  variables,  then  this  exclusion  does  not  change 
the  failure  probability.  Our  convention,  however,  facilitates  easy  transfer  of  the  results  in  [20]  to 
the  framework  of  the  present  paper  and  therefore  allows  general  forms  of  the  limit-state  function 
and  a  wide  range  of  probability  distributions.  Second,  while  we  define  g{x,v)  >  0  as  failure,  some 
authors  adopt  the  opposite  convention  where  g{x,  v)  <  0  represents  failure.  Obviously,  it  is  trivial 
to  switch  between  the  two  conventions  by  multiplying  the  limit-state  function  with  —1.  In  this 
paper,  we  use  the  convention  g{x,  v)  >  0  to  indicate  failure  as  our  derivations  appear  simpler  in 
that  case. 

In  Section  2,  we  discuss  the  properties  of  the  failure  probability  in  detail.  Section  3  presents 
the  buffered  failure  probability  and  shows  that  it  is  more  conservative  than  the  failure  probability, 
accounts  for  unlikely  but  possible  realizations  of  the  state  of  the  structure,  and  has  significant 
computational  advantages.  Section  4  generalizes  the  discussion  to  structural  systems  with  multiple 
limit-state  functions.  Section  5  illustrates  the  use  of  the  buffered  failure  probability  in  design 
optimization  of  a  truss  structure  and  a  vehicle  frame.  We  end  the  paper  with  concluding  remarks 
in  Section  6. 
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2  Properties  of  the  Failure  Probability 


While  the  definition  of  the  failure  probability  is  appealing  due  to  its  relative  simplicity,  it  exhibits 
several  undesirable  properties  resulting  in  significant  theoretical  and  practical  difficulties.  We  dis¬ 
cuss  these  in  turn  next. 

2.1  Simplistic  Characterization  of  Structures  as  Failed  or  Safe 

The  current  approach  to  structural  reliability  effectively  characterizes  a  structure  to  be  in  only 
one  of  two  possible  states:  failed,  i.e.,  5'(x,  v)  >  0,  or  safe,  i.e.,  g{x,'v)  <  0.  Consequently,  the 
“degree”  of  failure  is  not  important.  For  example,  the  event  {^'(x,  V)  =  100}  is  no  worse  than 
the  event  {g'(x,  V)  =  0.01}  as  they  both  are  subsets  of  the  failure  domain  and  contribute  to  the 
failure  probability.  However,  a  designer  would  most  likely  prefer  the  event  {^(x,  V)  =  0.01}  as  it 
represents  only  a  minor  violation  of  a  threshold,  possibly  somewhat  arbitrarily  set.  On  the  other 
hand,  the  event  {g'(x,  V)  =  100}  may  be  catastrophic.  The  theory  of  structural  reliability  does  not 
account  for  the  designer’s  preference  in  this  case.  This  preference  may  become  important  when  a 
designer  compares  two  candidate  designs  as  the  following  example  illustrates. 

Example  1.  Consider  the  design  of  a  structure  that  is  characterized  by  the  limit-state  function 

^(x,  v)  =  100  -  xiui  -  (1  -  xi)u2,  (2) 

where  100  is  a  deterministic  load  on  the  structure  and  xi  is  a  design  variable  to  be  chosen  by  the 
designer.  Only  xi  =  0  and  xi  =  1  are  allowable  choices.  Moreover,  let  Vi  be  a  normally  distributed 
random  variable  with  mean  150  and  standard  deviation  15  representing  the  strength  of  the  structure 
when  design  xi  =  1.  When  design  xi  =  0,  the  strength  of  the  structure  is  V2  which  is  a  random 
variable  with  mean  150  and  a  triangular  probability  density  function  in  the  range  [98.40,  175.8] 
with  values  near  175.8  being  the  most  likely  outcomes.  Figure  1  illustrates  the  probability  density 
functions  of  5^(0,  V)  and  gf(l,  V).  For  both  designs,  the  probability  of  failure  is  4.29- 10“^.  However, 
as  seen  from  Figure  2,  which  depicts  the  upper  tails  of  the  probability  density  functions  in  Figure 
1,  the  probability  of  an  “extreme  event”  is  substantial  in  case  of  design  xi  =  1,  but  nonexistent 
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Figure  1:  Example  1:  Probability  density  functions  (pdf)  of  ^((1,  V)  and  ^'(0,  V). 

for  design  xi  =  0.  For  example,  the  probability  of  the  event  {^(x,  V)  >  2}  is  2.63  •  10“^  for  design 
xi  =  1  but  for  design  xi  =  0  that  probability  is  of  course  zero.  While  this  is  obviously  an  artificial 
example,  it  illustrates  that  two  designs  with  the  same  failure  probability  may  have  significantly 
different  characteristics.  If  the  designer  only  computes  the  failure  probability,  this  difference  may 
not  be  revealed. 

2.2  Inaccurate  or  Computationally  Costly  Approximations 

Since  the  uncertainty  in  a  structure  often  needs  to  be  characterized  by  many  (hundreds  of)  random 
variables,  the  computation  of  the  failure  probability  for  a  given  design  x  requires  the  evaluation 
of  a  high-dimensional  integral,  see  (1).  As  that  evaluation  is  usually  impossible  to  carry  out 
analytically  and  computationally  expensive  to  carry  out  by  numerical  integration,  approximations 
based  on  Monte-Carlo  simulation  and  geometric  considerations  are  typically  used. 
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Figure  2:  Example  1;  Tails  of  probability  density  functions  (pdf)  of  ^((1,  V)  and  ^'(0,  V). 

2.2.1  Monte-Carlo  Simulation 

For  a  given  design  x,  Monte-Carlo  simulation  estimates  the  failure  probability  p(x)  by  generating 
N  independent  realizations  v^,  v^,  of  the  random  vector  V  and  computing  the  fraction  of 

the  realizations  in  the  failure  domain,  i.e.,  the  estimate  of  the  failure  probability 

1  ^ 

P7v(x)  =  —  ^/(5(x,vJ)  >  0).  (3) 

i=i 

The  corresponding  estimator  is  unbiased  and,  from  the  central  limit  theorem,  we  know  that  the 
standard  deviation  of  the  estimator  decays  proportional  to  1/y/N,  as  N  ^  oo.  While  this  decay 
rate  cannot  be  improved  upon,  the  standard  deviation  of  the  estimator  can  often  be  much  improved 
by  the  use  of  variance  reduction  techniques  such  as  importance  sampling,  directional  sampling,  and 
Markov-Chain  Monte-Carlo  sampling,  see,  e.g.,  [26].  Since  the  standard  deviation  decays  only 
proportional  to  1/y/N  and  the  effort  required  to  compute  pn{^)  grows  linearly  in  N ,  Monte  Carlo 
simulation  is  usually  computationally  costly. 
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2.2.2  Geometric  Approximations 


If  the  random  vector  V  consists  of  independent  standard  normal  random  variables  and  the  limit- 
state  function  is  affine  in  v,  i.e.,  g{x,v)  =  a(x)'v  -|-  6(x)  for  some  m-valued  function  a(x)  and 
real-valued  function  6(x),  then  the  failure  probability  p(x)  =  <I>(— /3(x))  whenever  p(x)  <  0.5, 
see,  e.g..  Chapters  4-6  in  [5].  Here,  <!>(•)  is  the  cumulative  distribution  function  of  a  standard 
normal  random  variable  and  /I(x)  is  the  shortest  distance  from  the  origin  in  IR”^  (he.,  the  space  of 
realizations  of  V)  to  the  surface  {v|gf(x,  v)  =  0},  see  Figure  3,  where  5^1  (x,v)  is  an  example  of  an 
affine  limit-state  function.  We  refer  to  /3(x)  as  the  reliability  index  of  design  x.  It  can  be  shown 
that  in  this  case 

/3(x)  = -6(x)/||a(x)||.  (4) 

If  5'(x,  v)  is  not  affine,  see,  e.g.,  5(2(x,  v)  in  Figure  3,  then  <I>(— /3(x))  is  an  approximation  of  the 
failure  probability.  In  this  case,  there  is  no  explicit  expression  for  /3(x)  and  it  must  be  computed 
by  solving  the  optimization  problem 

/3(x)  =  min  ||v||  (5) 

s.t.  g{x,  v)  >  0. 

There  is  empirical  evidence  that  the  approximation  $(— /3(x))  of  the  failure  probability  is  quite  ac¬ 
curate  on  classes  of  applications  arising  in  structural  engineering;  see  for  example  [36]  and  references 
therein.  However,  the  approach  may  also  lead  to  inaccuracy  as  discussed  below. 

When  5((x,  v)  is  not  concave^  in  v,  this  optimization  problem  may  have  points  satisfying 
the  Karush-Kuhn- Tucker  (KKT)  first-order  necessary  conditions  for  a  local  minimum  but  that 
are  not  global  minima.  For  example,  limit-state  function  gf3(x,  v)  in  Figure  3  results  in  a  line 
{v|(5r3(x,  v)  =  0}  with  many  points  that  are  locally,  but  not  globally,  the  closest  point  to  the  origin. 
Since  standard  nonlinear  optimization  algorithms  such  as  SNOPT  [9],  LANCELOT  [4],  and  NLPQL 
[28]  only  guarantee  convergence  to  such  a  KKT  point,  it  may  be  difficult  to  compute  the  globally 
optimal  solution  of  (5)  in  this  situation,  let  alone  prove  that  an  obtained  point  is  globally  optimal. 

The  same  holds  for  algorithms  specialized  for  solving  (5)  such  as  the  iHLRF  algorithm  [13].  Hence, 
^See,  e.g.,  [2,  3]  for  definitions  of  concavity  and  convexity. 
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Figure  3:  Reliability  indices  /3(x)  for  three  limit-state  functions  and  a  given  design  x.  The  shaded 
areas  indicate  failure  domains. 

/?(x)  could  be  significantly  overestimated,  thereby  leaving  serious  design  risks  undetected.  For 
example,  a  standard  nonlinear  programming  algorithm  may  return  the  same  value  for  the  three 
limit-state  functions  in  Figure  3  when  applied  to  (5).  The  value  would  be  correct  for  gfi(x,  v)  and 
g'2(x,  v),  but  severely  overestimate  the  reliability  index  for  g'3(x,  v).  Even  if  the  global  minimum  is 
found  in  (5),  we  see  from  Figure  3  that  d>(— /3(x))  may  overestimate  p(x),  as  in  the  case  of  52(x,  v), 
or  underestimate  it  as  in  the  case  of  v).  In  general,  it  is  difficult  to  know  how  close  <h(— /3(x)) 
is  to  p(x). 

In  practice,  V  is  essentially  never  a  vector  of  independent  standard  normal  random  variables. 
Hence,  to  apply  the  above  approximation  one  typically  needs  to  carry  out  a  probability  transforma¬ 
tion,  see,  e.g..  Chapter  7  of  [5].  Random  vectors  governed  by  distributions  such  as  the  multivariate 
normal  (possibly  with  correlation)  and  lognormal  distributions  can  be  transformed  into  a  standard 
normal  vector  using  a  smooth  bijective  mapping.  Other  transformations  can  also  be  carried  out  at 
least  approximately.  A  transformation  can  make  the  limit-state  function  highly  nonlinear  and  non- 
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concave  as  function  of  the  independent  standard  normal  random  variables,  which  makes  it  difficult 
to  determine  the  global  minimum  of  (5). 

The  method  of  estimating  p(x)  by  <h(— /?(x))  is  referred  to  as  the  hrst-order  reliability  method 
as  it  effectively  linearizes  a  transformed  limit-state  function.  An  extension  of  this  method  is  the 
second-order  reliability  method  where  the  transformed  limit-state  function  is  approximated  by  a 
quadratic  function,  see  Chapter  6  of  [5].  However,  the  second-order  reliability  method  suffers  from 
the  same  difficulties  as  the  first-order  method,  though  its  accuracy  may  be  better.  An  alternative 
method  is  to  attempt,  after  a  transformation  to  independent  standard  normal  random  variables, 
to  determine  the  largest  ball  in  IR™,  centered  at  the  origin,  with  g'(x,  v)  <  0  for  all  v  in  the  ball. 
Using  the  chi-square  distribution,  this  leads  to  an  upper  bound  on  the  failure  probability  p(x). 
However,  the  bound  is  usually  overly  conservative  and  of  little  practical  use. 

2.3  Poorly  Behaving  Sensitivities  of  Failure  Probability  and  its  Approximations 

In  sensitivity  analysis  and  design  optimization,  we  examine  the  effect  on  the  failure  probability 
(or  its  approximation)  of  inhnitesimal  changes  in  the  design.  Hence,  differentiability  of  the  failure 
probability  and  its  approximations  with  respect  to  design  x  as  well  as  computable  formulae  for  the 
corresponding  gradient  become  important.  Specifically,  standard  nonlinear  optimization  algorithms 
require  all  functions  in  an  optimization  problem  to  be  continuously  differentiable.  If  this  condition 
is  not  satisfied,  the  algorithms  may  break  down  without  returning  an  optimized  design. 

2.3.1  Gradient  of  the  Failure  Probability 

The  issue  of  differentiability  of  the  failure  probability  is  nontrivial  as  the  integrand  in  (1)  is  not 
differentiable.  (The  indicator  function  makes  a  jump  from  1  to  0  as  the  condition  g{x,  v)  >  0  goes 
from  being  satisfied  to  not  satisfied.)  Hence,  we  cannot  simply  compute  the  derivative  of  an  integral 
by  integrating  the  derivative  of  the  integrand  which  is  allowed  under  weak  assumptions  when  the 
integrand  is  differentiable. 

Despite  this  situation,  the  failure  probability  is  actually  continuously  differentiable  with  respect 
to  the  design  x  under  rather  general  conditions  when  the  failure  domain  is  bounded  and  the  limit- 
state  function  is  continuously  differentiable  with  respect  to  the  design  [33].  However,  the  gradient 
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formula  in  [33]  is  difficult  to  use  in  estimation  because  it  may  involve  surface  integrals.  In  [14] 
(see  also  [15]),  an  integral  transformation  is  presented,  which,  when  it  exists,  leads  to  a  simple 
formula  for  the  gradient  of  the  failure  probability.  However,  it  is  not  clear  under  what  conditions 
the  transformation  exists.  As  in  [33],  [32]  assumes  that  the  failure  domain  is  bounded.  With  this 
restriction  as  well  as  the  assumption  that  the  failure  domain  is  “star-shaped,”  a  formula  for  the 
gradient  of  the  failure  probability  involving  integration  over  a  simplex  is  derived.  In  principle,  this 
integral  can  be  evaluated  by  Monte  Carlo  simulation.  However,  to  the  authors’  knowledge,  there  is 
no  computational  experience  with  estimation  of  failure  probabilities  for  highly  reliable  mechanical 
structures  using  this  formula. 

In  Section  9.2  of  [5],  with  generalizations  and  proofs  in  [24],  we  find  convenient  expressions  for 
the  gradient  of  the  failure  probability  under  similar  assumptions  to  those  in  [32].  The  expressions 
can  be  estimated  using  Monte  Carlo  simulation  with  good  accuracy  at  moderate  computational  ex¬ 
pense  when  the  star-shaped  assumption  is  satisfied  and  the  number  of  random  variables  is  moderate. 
However,  it  becomes  increasingly  costly  to  estimate  the  expression  using  Monte  Carlo  simulation 
when  the  number  of  random  variables  grows.  Moreover,  in  practice,  it  is  difficult  to  verify  the 
star-shape  assumption.  An  alternative  formula  for  the  gradient  of  the  failure  probability  is  pre¬ 
sented  in  [23,  22]  that  can  also  be  estimated  using  Monte  Carlo  simulation.  However,  the  formula 
relies  on  the  implicit  function  theorem  applied  to  the  equation  g{x,  v)  =  0  that  may  not  always  be 
applicable. 

2.3.2  Gradient  of  the  Reliability  Index 

As  described  in  Subsection  2.2,  the  failure  probability  p(x)  can  rarely  be  computed  exactly  and  the 
approximation  <I>(— /3(x))  is  often  used,  where  the  reliability  index  /3(x)  is  defined  in  (5).  Since  the 
cumulative  distribution  function  <!>(•)  is  continuously  differentiable,  differentiability  of  this  approx¬ 
imation  depends  on  the  properties  of  /3(x).  We  find  expressions  for  the  gradient  of  /3(x)  in  Chapter 
8  of  [5],  but  those  cannot  hold  for  all  x  as  the  following  simple  example  illustrates. 
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Figure  4:  Example  2.  Equation  gf(x,  v)  =  v\lx\  +  ul  —  1  =  0  illustrated  for  three  different  values 
of  xi-  The  shaded  areas  indicate  the  failure  domain. 

Example  2.  Consider  the  limit-state  function 

v‘i 

5(x,v)  =  ^ -  1,  (6) 

let  Vi  and  V2  be  independent  standard  normal  random  variables,  and  let  xi  >  0  be  a  design  variable, 
see  Eigure  4.  As  /?(x)  is  defined  as  the  distance  to  the  closest  point  on  the  surface  {v|5r(x,  v)  =  0}, 
see  (5),  we  find  that  /3(x)  =  xi  if  0  <  xi  <  1  and  1  otherwise.  Hence,  dl3{x.)/dxi  =  1  if  0  <  xi  <  1, 
dl3{x)/dxi  =  0  if  xi  >  1,  and  the  derivative  is  not  defined  when  xi  =  1,  see  Eigure  5.  As  we  see 
from  this  figure,  /?(x)  is  not  continuously  differentiable  and  the  derivative  at  xi  =  1  is  not  defined. 

As  Example  2  illustrates,  /3(x)  may  not  be  continuously  differentiable  and,  hence,  standard 
nonlinear  optimization  algorithms  may  stall  at  points  that  are  not  KKT  points  when  applied  to 
design  optimization  models  involving  /3(x). 

In  view  of  the  above  discussion,  we  see  that  the  differentiability  of  the  failure  probability  as 
well  as  the  existence  of  tractable  formulae  for  its  gradient  rely  on  assumptions  that  may  not  hold 
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Figure  5;  Example  2.  Illustration  of  reliability  index  /3(x)  (solid  line)  and  its  derivative  (dashed 
line). 

and  that  are  difficult  to  verify  in  practice.  Moreover,  the  frequently  used  reliability  index  provides 
an  approximation  of  the  failure  probability  <h(— /3(x))  that  may  not  be  continuously  differentiable. 
Hence,  even  if  the  limit-state  function  is  a  continuously  differentiable  function  in  the  design  vari¬ 
ables,  the  failure  probability  and  <!>(— /3(x))  may  not  be. 

2.4  Lack  of  Convexity  of  the  Failure  Probability 

As  stated  above,  standard  nonlinear  optimization  algorithms  typically  only  guarantee  convergence 
to  a  KKT  point.  However,  if  a  design  optimization  problem  has  a  convex  objective  function,  which 
we  would  like  to  minimize,  and  the  constraints  form  a  convex  feasible  region,  then  a  KKT  point 
must  be  a  global  optimal  design  for  the  problem.  Absent  convexity,  it  may  be  difficult  to  compute 
a  globally  optimal  design,  let  alone  prove  that  an  obtained  design  is  globally  optimal.  Therefore 
we  would  like  to  formulate  convex  design  optimization  models  if  possible.  We  refer  to  [3]  for  an 
introduction  to  convex  optimization. 

From  this  discussion  we  conclude  that  the  convexity  of  the  failure  probability  p(x)  would  be 
valuable  when  solving  a  design  optimization  problem.  Unfortunately,  it  is  unknown  whether  p(x) 
is  convex  even  if  g{x,  v)  is  convex  in  x  for  all  v.  Hence,  p(x)  does  not  “preserve”  convexity  as  a 
convex  limit-state  function  may  result  in  a  nonconvex  failure  probability.  The  same  situation  holds 
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when  the  reliability  index  /3(x)  is  used  to  approximate  the  failure  probability.  For  this  reason, 
we  expect  that  design  optimization  problems  involving  the  failure  probability  or  the  reliability 
index  may  have  many  local  minima  that  are  not  globally  optimal.  Standard  nonlinear  optimization 
algorithms  are  unlikely  to  find  the  globally  optimal  design  and  may  return,  at  best,  locally  optimal 
designs.  Consequently,  it  may  be  necessary  to  apply  computationally  expensive  global  optimization 
algorithms,  see,  e.g.,  [10]. 

3  Buffered  Failure  Probability 

As  reviewed  in  Section  2,  the  failure  probability  has  several  troublesome  properties.  In  this  section, 
we  discuss  an  alternative  probability,  which  we  call  the  buffered  failure  probability,  that  has  several 
advantages  over  the  failure  probability.  The  buffered  failure  probability  relates  to  the  conditional 
value-at-risk  [19,  20],  which  is  now  widely  used  in  the  area  of  financial  engineering  to  assess  invest¬ 
ment  portfolios.  The  tutorial  paper  [18]  provides  an  overview  including  relation  to  safety  margins 
and  potential  replacements  for  failure  probability  constraints.  However,  buffered  failure  probability 
is  directly  introduced  and  explained  here  for  the  first  time. 

3.1  Definition 

We  first  recall  that  for  any  probability  level  a,  the  a-quantile  of  the  distribution  of  a  random 
variable  is  the  value  of  the  inverse  of  the  corresponding  cumulative  distribution  function  at  a. 
For  simplicity  in  presentation,  we  assume  here  and  throughout  this  paper  that  the  cumulative 
distribution  function  of  g{x,  V)  is  continuous  and  strictly  increasing  for  all  x.  For  definitions  which 
serves  to  fully  generalize  beyond  this  case,  we  refer  to  [20].  We  consider  especially  the  random 
variable  (/(x,  V)  for  a  given  design  x  and  denote  the  a-quantile  of  5'(x,  V)  by  (?a(x).  As  indicated 
by  the  notation,  ^^(x)  depends  on  the  design  x  as  the  probability  distribution  of  gr(x,  V)  changes 
with  X.  Figures  6  and  7  illustrate  ^a(x)  for  the  case  when  g{x,\')  is  normally  distributed  with 
mean  —1  and  standard  deviation  1.  Figure  6  shows  the  cumulative  distribution  function  of  gf(x,  V) 
in  this  case  and  quantiles  corresponding  to  probability  levels  a  =  0.60  and  ag  =  0.84.  Figure 
7  illustrates  the  same  information  using  the  probability  density  function  of  g(x,  V)  and,  hence. 
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Figure  6:  Cumulative  distribution  function  (cdf)  of  with  examples  of  a-quantile  (?a(x)  and 

a-superquantile  5q(x)  when  normally  distributed  with  mean  —1  and  standard  deviation  1. 
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Figure  7:  Probability  density  function  (pdf)  of  g{x,Y)  with  examples  of  a-quantile  ^^(x)  and 
a-superquantile  5q(x)  when  normally  distributed  with  mean  —1  and  standard  deviation  1. 
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probabilities  correspond  to  areas  under  that  function.  In  view  of  Figures  6  and  7  and  (1),  we 
find  that  the  failure  probability  is  equal  to  one  minus  the  probability  level  that  results  in  the 
quantile  being  zero.  For  example,  in  Figure  6  we  find  that  ao  =  0.84  gives  qaoi^)  =  0.  Hence, 
p{x)  =  l-ao  =  l-  0.84  =  0.16. 

Before  we  dehne  the  buffered  failure  probability,  we  introduce  a  quantity  that  is  closely  related 
to  the  quantile.  For  any  probability  level  a,  we  define  the  a-superquantile  as 


qa{x)  =  E[g{x,V)\g{x,Y)  >  ga(x)], 


(7) 


where  the  vertical  bar  indicates  a  conditional  expectation.  That  is,  the  a-superquantile  is  the 
average  value  of  g{x,Y),  conditional  on  the  event  that  g{x,Y)  is  no  less  than  the  a-quantile.  This 
quantity  is  called  Conditional  Value-at-Risk  in  financial  engineering,  but  we  here  propose  and  adopt 
the  application-independent  name  superquantile.  Figures  6  and  7  illustrate  the  superquantiles  of 
g{x,Y)  for  probability  levels  a  =  0.60  and  ao  =  0.84.  Since  g{x,Y)  is  normally  distributed,  it  is 
trivial  to  compute  the  superquantiles  using  the  well-known  conditional  expectation  formula  (see. 


e.g.,  [35]), 


9a  (x) 


P  + 


0'</>(ga) 

1  —  a  ’ 


(8) 


for  a  normally  distributed  (/(x,  V)  with  mean  g,  standard  deviation  a,  and  truncation  level  qa, 
where  (/>(•)  is  the  standard  normal  probability  density  function  and  qa  is  the  a-quantile  of  the 
standard  normal  distribution.  When  g{x,  V)  is  not  normally  distributed,  the  calculation  of  the 
superquantile  appears  much  more  difficult.  As  seen  in  the  next  subsection,  however,  it  can  be 
computed  in  a  remarkably  efficient  manner. 

Figure  7  highlights  the  dehnition  of  a  superquantile  as  a  conditional  expectation.  As  seen  for 
probability  level  a  =  0.60,  the  corresponding  quantile  is  —0.75.  The  corresponding  superquantile 
is,  roughly  speaking,  the  value  that  splits  the  interval  [— 0.75,oo)  into  two  “balancing”  parts.  The 
area  under  the  probability  density  function  between  —0.75  and  the  value  (the  lightly  shaded  area  in 
Figure  7)  “balances”  the  area  under  the  function  above  the  value  (the  heavily  shaded  area).  In  this 
case,  that  value  is  —0.03  as  computed  by  (8).  Similarly,  for  probability  level  0.84,  the  corresponding 
quantile  is  0  and  the  corresponding  superquantile  is  0.53.  That  is,  the  area  under  the  probability 
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density  function  between  0  and  0.53  “balances”  the  area  under  the  function  above  0.53. 


We  note  that  in  general  (7a (x)  <  qa{^)  for  any  probability  level  a  and  design  x.  In  [18]  we  also 
find  the  following  equivalent  formula  for  the  superquantile: 

1 


9a  (x)  = 


1  —  a  . 


qa'{^)da' . 


(9) 


We  do  not  repeat  the  derivation  of  this  expression  here,  but  note  that  the  expression  essentially 
averages  the  quantiles  for  probability  levels  larger  than  a. 

We  now  define  the  buffered  failure  probability  p(x)  to  be  equal  to  1  —  a  where  a  is  selected 
such  that  the  superquantile 

ga(x)  =  0.  (10) 


That  is, 

p(x)  =  P[ff(x,  V)  >  ga(x)],  (11) 

where  a  is  selected  such  that  (10)  holds.  Hence,  9i-p(x)(x)  =  0.  We  see  from  Figures  6  and  7 
that  the  probability  levels  a  =  0.60,  which  led  to  ga(x)  =  —0.03,  and  ao  =  0.84,  which  led  to 
9qo(x)  =  0.53,  are  slightly  too  small  and  much  to  large,  respectively,  to  result  in  a  corresponding 
superquantile  of  zero.  However,  it  is  easy  to  find  by  trial-and-error  and  (8)  that  a  probability  level 
a  =  0.62  results  in  a  quantile  of  —0.70  and  a  superquantile  of  approximately  zero  as  illustrated 
in  Figures  8  and  9.  (We  present  a  much  easier  way  than  trial-and-error  below  for  computing  the 
superquantile.)  By  definition,  see  (11),  the  buffer  probability  is  then  1  —  a  =  1  —  0.62  =  0.38,  which 
is  somewhat  larger  than  the  failure  probability  of  0.16. 

In  general,  we  find  that 

p(x)  <  p(x)  (12) 

for  any  x,  see  [19,  20,  18].  Hence,  the  buffered  failure  probability  is  a  conservative  estimate  of  the 
failure  probability  for  any  design  x.  As  we  see  below,  the  degree  of  overestimation  is  usually  modest. 
We  stress,  however,  that  the  buffered  failure  probability  carries  more  information  about  the  design 
than  the  failure  probability  as  it  includes  information  about  the  upper  tail  of  (^(x,  V).  Hence, 
for  designs  where  the  probability  of  (7(x,  V)  taking  on  values  substantially  above  zero  is  relatively 
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Figure  8:  Cumulative  distribution  function  (cdf)  of  gr(x,  V),  as  in  Figure  6,  with  a  selected  such 
that  the  a-superquantile  ga(x)  =  0.  Illustration  of  the  buffered  failure  probability  p(x)  and  the 
failure  probability  p(x). 

large,  the  buffered  failure  probability  tends  to  be  somewhat  larger  than  the  failure  probability. 
In  contrast,  if  the  probability  of  g(x,  V)  taking  on  large  values  is  small,  then  the  buffered  failure 
probability  is  typically  close  to  the  failure  probability. 

As  we  discuss  below,  the  buffered  failure  probability  is  surprisingly  easy  to  compute,  possesses 
several  convenient  properties,  and  avoids  many  of  the  difficulties  associated  with  the  failure  prob¬ 
ability.  Hence,  we  believe  there  are  substantial  advantages  to  replacing  the  failure  probability  by 
the  buffered  failure  probability  in  engineering  design. 

Example  3.  Consider  the  limit-state  function  given  in  Example  1  and  recall  that  p(0)  =  p{l)  = 
4.29  •  10“^.  We  now  compute  the  buffered  failure  probability  for  the  designs  xi  =  0  and  xi  =  1. 
Since  5(0,V)  is  given  by  a  triangular  probability  density  function,  we  determine  an  ol  such  that 
^o(O)  =  0  by  integration  and  find  that  p(0)  =  1  —  a  =  9.65  •  10“^.  For  design  xi  =  1,  go(x) 
is  the  expectation  of  a  truncated  normal  distribution,  which  is  easily  calculated  by  (8).  We  use 
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Figure  9;  Probability  density  function  (pdf)  of  gr(x,  V),  as  in  Figure  7,  with  a  selected  such  that 
the  a-superquantile  gQ(x)  =  0.  Illustration  of  the  buffered  failure  probability  p(x)  and  the  failure 
probability  p(x). 
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trial-and-error  to  determine  an  a  such  that  5a(l)  =  0  and  find  that  p{l)  =  1  —  a  =  1.13  •  10“^. 
We  first  observe  that  both  designs  satisfy  (12)  as  expected.  We  also  see  that  design  xi  =  0  has  a 
smaller  buffered  failure  probability  than  design  xi  =  1  and  is  therefore  “safer”  in  the  sense  of  the 
buffered  failure  probability.  This  corresponds  to  our  intuition  discussed  in  Example  1,  where  we 
concluded  that  design  xi  =  0  was  preferable  due  to  the  smaller  probability  of  extreme  violation  of 
the  threshold. 

Example  3  illustrates  the  fact  that  the  buffered  failure  probability  takes  into  account  the  tail 
behavior  of  the  distribution  of  5'(x,  V)  and  hence  offers  an  alternative  measure  of  reliability  of  a 
structure  that  may  better  reflect  designers’  concerns. 

As  seen  from  (11)  and  (1),  the  buffered  failure  probability  shifts  the  threshold  level  from  zero 
downwards  to  qai^)  (a  negative  number)  and  therefore  adds  a  “buffer  zone”  to  the  failure  domain. 
We  observe  that  the  threshold  shift  and  the  buffer  zone  depend  on  the  probability  distribution  of 
g{x,  V)  and,  hence,  on  x  as  illustrated  in  Eigure  10  for  the  limit-state  function  in  Examples  1  and  3. 
That  figure  shows  two  solid  and  two  dotted  lines.  The  vertical  solid  line  represents  {v  |  gf(l,  v)  =  0}, 
i.e.,  the  limit  between  the  failure  domain  (to  the  left)  and  its  complement  the  safe  domain  (to  the 
right)  for  design  xi  =  1.  The  dotted  vertical  line  represents  {v  |  g{l,v)  =  ^^(l)}.  This  line 
shifts  to  the  right  as  compared  to  the  solid  line  resulting  in  a  buffer  zone  (shaded  dark).  Similarly, 
the  horizontal  solid  line  represents  {v  |  g(0,v)  =  0},  i.e.,  the  limit  between  the  failure  domain 
(below)  and  the  safe  domain  (above)  for  design  xi  =  0.  The  dotted  horizontal  line  represents 
{v  I  gf(0,  v)  =  (7o(0)}.  This  line  shifts  up  as  compared  to  the  solid  line  resulting  in  a  narrow  buffer 
zone.  We  observe  that  the  threshold  shift  and  buffer  zone  are  substantially  smaller  for  xi  =  0  than 
for  xi  =  1.  In  general,  the  line,  surface,  or  hypersurface  {v  |  5f(x,  v)  =  ga(x)}  may  not  be  parallel 
to  {v  I  5r(x,  v)  =  0}  for  a  given  x.  In  Eigure  10,  however,  those  lines  are  parallel  due  to  the  fact 
that  the  limit-state  function  in  Examples  1  and  3  is  affine  in  v. 

Erom  the  above  definition  of  the  superquantile,  it  may  appear  difficult  to  compute  the  buffered 
failure  probability  in  general.  However,  this  is  not  the  case  as  the  next  subsection  describes. 
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Figure  10:  Contours  of  limit-state  function  in  Examples  1  and  3  for  xi  =  0  and  1. 

3.2  Buffered  Failure  Probability  in  Design  Optimization 

Suppose  we  would  like  find  a  design  with  failure  probability  no  larger  than  a  threshold  1  —  uq-  That 
is,  we  would  like  to  determine  a  design  x  that  satisfies  the  constraint 

p(x)  <  1  -  ao-  (13) 


In  view  of  Section  2,  we  observe  that  standard  optimization  algorithms  may  have  substantial 
difficulties  on  problems  with  constraints  of  the  form  (13).  We  now  show  that  the  alternative 
constraint 

p(x)  <  1  -  ao  (14) 

in  terms  of  the  buffered  failure  probability  is  much  easier  to  handle.  We  start  by  noting  that  a 
design  x  that  satisfies  (14)  also  satisfies  (13).  Hence,  (14)  is  a  conservative  requirement. 

The  ease  with  which  (14)  can  be  handled  in  optimization  algorithms  clearly  hinges  on  our 
ability  to  evaluate  p(x)  or  equivalent  expressions.  While  p(x)  cannot  be  expressed  explicitly,  there 
is  a  convenient,  equivalent  expression  for  (14)  that  we  derive  next. 
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In  view  of  Figure  8,  we  see  that  (14)  holds  if  and  only  if 

gao(x)<0.  (15) 

It  is  shown  in  [19]  that 

g„(x)  =  minr/Q,(2;o,x),  (16) 

where  zq  is  an  auxiliary  design  variable  and 

r/Q,(zo,x)  =  zo  +  — — F;[max{0,5((x,  V)  -  zq}].  (17) 

I  —  a 

We  do  not  include  a  derivation  of  this  expression  as  it  is  somewhat  involved  and  refer  the  interested 
reader  to  [19].  Hence,  the  task  of  finding  a  design  x  that  satisfies  p(x)  <  1  —  ao  is  equivalent  of 
finding  a  design  x  and  an  auxiliary  variable  zq  such  that 

r7„o(zo,x)  <  0.  (18) 

Suppose  that  the  goal  is  to  determine  a  design  x  that  minimizes  some  continuously  differ¬ 
entiable  objective  function  /(x)  (e.g.,  cost)  subject  to  the  reliability  constraint  p(x)  <  1  —  ao 
and  a  finite  number  of  continuously  differentiable  equality  and  inequality  constraints  abstractly 
represented  by  the  set  X.  That  is,  we  would  like  to  solve  the  design  optimization  problem 

P  :  min/(x) 

s.t.  p(x)  <  1  —  ao 

X  G  X. 


In  view  of  the  discussion  above,  the  alternative  formulation  in  terms  of  the  buffered  failure  proba¬ 
bility  takes  the  form 


BP  :  min  /(x) 

X,20 

S.t.  Zq - iii[max{0, 5((x,  V)  —  zo}]  <  0 

1  -  ao 

X  G  X, 


where  we  observe  that  the  optimization  is  over  both  x  and  zq. 
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We  usually  cannot  compute  £^[max{0,  ^((x,  V)  —  zq}]  explicitly.  However,  the  expectation  can 
be  estimated  by  its  sample  average.  Let  v^,...,v'^  be  realizations  of  V.  Then,  the  optimization 
problem 


BP' 


N 


min  /(x) 


X,ZO 

N 


S.t.  Zq  + 


^(1  -  ao) 


^max{0,5r(x,  v^)  -  zo}  <  0 


(19) 


X  G  X, 


is  an  approximation  of  the  problem  BP.  Even  if  the  limit-state  function  gf(x,  v)  is  continuously 
differentiable  for  all  v,  BP'^r  is  not  directly  tractable  by  standard  nonlinear  optimization  algorithms 
due  to  the  nonsmoothness  of  the  max-function  in  (19).  BP'^r  is  solvable  by  a  specialized  algorithm 
found  in  [16],  but  we  do  not  describe  that  algorithm  here.  Instead  we  present  an  equivalent 
transcription  of  BP'^r  that  facilitates  the  use  of  standard  nonlinear  optimization  algorithms. 

We  let  zi,  ...,Z]\f  be  auxiliary  design  variables  and  denote  z  =  {zo,zi,  ...,zm)' ■  Then,  BP'jv  is 
equivalent  to  the  following  intermediate  problem 


s.t.  zo  -I- 


min  /(x) 

X.Z 

1 

N 

0 

VI 

N{1  -  ao) 

:{0,c/(x,  v^) 

-  2:0} 

= 

X 

(20) 


where  we  simply  force  the  auxiliary  design  variables  to  take  on  the  “right”  values.  We  can  relax 
the  equality  constraints  to  less-than-or-equal  constraints  as  there  is  no  benefit  to  let  the  variables 
take  on  values  such  as  max{0, ^((x, v-^)  —  zq}  <  Zj  for  any  j  =  1,2,...,N.  Moreover,  a  constraint 
of  the  form  max{0, 5'(x,  v-^)  —  zq}  <  Zj  is  equivalent  to  the  two  constraints  g(x,v^)  —  zq  <  Zj  and 
0  <  Zj.  This  leads  to  the  following  equivalent  problem  of  BP'^y: 


BP/v  : 


min  /(x) 

X,Z 


s.t.  Zo  -I- 


1  ^ 

A'(l-«o)§ 


<  0 
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g{yL,w^)-zo  <  Zj,  j  =  l,2,...,N 
Zj  >  0,  j  =  1,2,  ...,7V 


(21) 


X  G  X. 

We  propose  that  engineers  consider  BP  tv  instead  of  P  when  designing  structures  for  reasons  sum¬ 
marized  next. 

3.3  Comparison  of  Probabilities 

There  are  four  main  advantages  to  consider  BPjv  instead  of  P.  First,  as  discussed  above,  the  failure 
probability  p{x)  and  its  gradient  cannot  generally  be  computed  exactly  and  must  be  approximated 
in  ways  which,  in  some  cases,  might  turn  a  blind  eye  to  serious  risks.  The  first-order  approxima¬ 
tion  <h(— /3(x))  has  unknown  accuracy  and  may  not  be  continuously  differentiable.  Monte  Carlo 
estimates  of  p(x)  have  error  bounds,  but  the  estimates  have  gradients  only  under  assumptions 
that  are  difficult  to  verify  in  practice.  Hence,  it  is  highly  problematic  to  apply  standard  nonlinear 
optimization  algorithms  to  optimization  problems  involving  p(x).  In  contrast,  BPtv  is  solvable  by 
standard  nonlinear  optimization  algorithms  as  long  as  the  limit-state  function  g{x,  v)  is  continu¬ 
ously  differentiable  with  respect  to  x.  This  is  a  substantially  less  stringent  condition  than  those 
required  for  P.  The  optimal  value  of  BPtv  is  close  to  the  optimal  value  of  BP  when  N  is  large 
(see  Chapter  4  of  [29]  for  specific  results  on  the  “proximity”  of  BPtv  to  BP).  Moreover,  BP  is  a 
restricted  problem  compared  to  P  because  the  buffered  failure  probability  overestimates  the  failure 
probability,  see  (12).  Hence,  a  feasible  design  in  BP  is  also  feasible  in  P. 

Second,  the  buffered  failure  probability  provides  an  alternative  measure  of  structural  reliability 
which  accounts  for  the  tail  behavior  of  the  distribution  of  (^(x,  V).  Hence,  designs  obtained  from 
BPtv  may  be  more  desirable  than  those  from  P. 

Third,  even  if  g{x,v)  is  convex  in  x,  p{x)  and  <h(— /3(x))  may  not  be  and,  hence,  it  may  be 
difficult  to  obtain  a  globally  optimal  design  of  P.  In  contrast,  the  region  defined  by  the  constraints 
(21)  is  convex  when  g{x,  v-^),  j  =  1,  2, ...,  N,  are  convex  functions  in  x.  Hence,  every  KKT  point  of 
BPtv  is  a  globally  optimal  design  when  /(x)  and  g{x,  v^),  j  =  1,2,  are  convex  functions  and 

the  region  X  is  a  convex  set.  Hence,  BPtv  “preserves”  convexity.  Even  if  not  all  of  these  conditions 
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are  satisfied,  we  expect  it  to  often  be  easier  to  determine  a  design  with  a  low  objective  function 
value  in  BP^r  than  in  P  because  BP^r  deals  with  g{x,v)  directly  instead  of  the  more  complex 
expression  p(x). 

We  expect  g{x,  v)  to  be  convex  in  x  in  several  practical  situations.  For  example,  suppose  that 
X  =  xi  represents  the  size  of  a  part  of  the  structure  and  the  strength  i?(x,  v)  of  the  structure 
grows  as  xi  grows  for  all  possible  realizations  v.  Moreover,  suppose  that  this  growth  in  strength 
is  constant  or  tapers  off  as  xi  grows.  Then,  i?(x,  v)  is  concave  for  all  v  and,  hence,  the  limit-state 
function  gr(x,  v)  =  ^(v)  —  i?(x,  v)  is  convex,  where  5(v)  describes  the  load  on  the  structure.  Since 
the  convexity  of  g{x,v)  with  respect  to  x  was  of  little  importance  in  the  context  of  p(x),  few  re¬ 
searchers  have  focused  on  developing  convex  limit-state  functions  or  approximations  thereof.  As  the 
importance  of  convexity  is  now  clear,  we  hope  that  this  paper  will  spur  research  into  the  develop¬ 
ment  of  convex  limit-state  functions.  While  physics  dictate  to  a  large  extent  the  form  of  limit-state 
functions,  engineers  may  still  have  opportunities  for  skillful  modeling,  including  the  development 
of  useful  approximations.  In  the  same  manner  as  a  simple  limit-state  function  g{x,v)  =  V2  —  xivi, 
which  is  linear  in  x,  is  equivalent  to  the  limit-state  function  g(x,v)  =  V2/{xivi)  —  1,  which  is  non¬ 
linear  in  X,  we  expect  the  development  of  (approximately)  equivalent  convex  limit-state  functions 
to  existing  nonconvex  limit-state  functions. 

Fourth,  BPtv  facilitates  the  development  of  approximation  schemes  for  limit-state  functions 
that  are  expensive  to  evaluate.  For  example,  if  the  evaluation  of  the  limit-state  function  involves 
the  output  of  a  finite  element  model,  it  may  not  be  possible  to  evaluate  the  limit-state  function 
more  than  a  few  hundred  or  a  few  thousand  times.  In  such  situations,  the  failure  probability  in  P 
is  often  replaced  by  response  surface  and  surrogate  models,  see,  e.g.,  [7,  31,  34].  This  allows  quick 
optimization,  but  the  quality  of  the  resulting  design  depends  on  the  fidelity  of  the  response  surface 
or  surrogate  model  used.  As  p(x)  may  be  a  highly  nonlinear,  nonconvex  function,  we  conjecture 
that  it  may  be  more  difficult  and  computationally  expensive  to  develop  a  good  surrogate  model  of 
p(x)  than  of  g{x,v),  about  which  we  may  have  problem-specihc  insight.  With  a  surrogate  model 
of  5((x,  v),  the  optimization  of  BPtv  using  that  surrogate  model  in  place  of  g{x,v)  can  often  be 
accomplished  relatively  quickly;  see  Section  5.2.  For  example,  suppose  that  x^,  k  =  l,2,...,iF,  is  a 
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selection  of  designs.  Then,  for  any  v  and  k, 


ff(x,  v)  PS  5(x^  v)  +  Vx5(x^  V)'(x  -  x'^)  (22) 

when  X  is  close  to  x^  and  g{'x.,  v)  is  continuously  differentiable  with  respect  to  x.  Obviously,  the 
selection  of  x^,  k  =  1,2, ...,  K,  e.g.,  by  means  of  an  experimental  design,  influences  the  accuracy  of 
this  approximation  and  is  an  important  topic  in  its  own  right.  In  this  paper,  however,  we  do  not 
discuss  this  topic  further.  Interested  readers  are  referred  to  [34]  and  references  therein.  Using  this 
linear  approximation  of  the  limit-state  function,  we  obtain  the  following  approximation  of  BPtv, 
which  is  intended  for  the  case  when  g{x,  v^)  is  convex  in  x  for  all  j  =  1,  2, ...,  N: 


LBPjv  :  min/(x) 

X,Z 


1  ^ 

< 

0 

+  Vx5'(x'',v^)'(x-x'')  -  Zo 

< 

Zj, 

j  =  l,2,...,iV,fc  =  l,2 

U' 

> 

0, 

j  =  l,2,...,N 

X 

G 

X. 

Under  that  convexity  assumption,  LBPjv  can  be  made  to  approximate  BP^r  arbitrarily  well  by 
selecting  more  designs  appropriately,  i.e.,  increasing  K.  We  note,  however,  that  LBPjv  is  a  noncon¬ 
servative  approximation  of  BP^r.  The  construction  of  conservative  approximations  of  BPjv  would 
also  be  possible  under  suitable  assumptions,  but  that  topic  is  beyond  the  scope  of  the  current  paper. 

Solving  LBPjv  only  requires  the  evaluation  of  the  limit-state  function  and  its  gradient  KN 
times  to  generate  the  problem  data  in  LBPjy.  During  optimization  no  evaluation  of  the  limit- 
state  function  or  its  gradient  is  needed  and,  hence,  can  be  carried  out  quickly.  If  the  objective 
function  /(x)  and  the  constraints  defining  X  are  linear,  then  LBPtv  is  a  linear  program  that  can 
be  solved  quickly  by  standard  linear  programming  solvers  or  decomposition  algorithms.  In  this 
case,  the  introduction  of  integrality  restrictions  on  x  may  also  be  tractable  as  this  makes  LBP a 
mixed-integer  linear  optimization  problem  that  often  can  be  solved  in  moderate  computing  times. 
In  comparison,  it  is  difficulty  to  solve  P  in  the  case  of  integrality  constraints  as  it  then  becomes  a 
mixed-integer,  nonlinear,  nonconvex,  optimization  problem. 
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3.4  Variance  Reduction 


While  a  large  sample  size  N  provides  a  good  approximation  of  BP^r  to  BP,  the  number  of  con¬ 
straints  and  decision  variables  in  BP  tv  grows  linearly  in  N.  The  accuracy  of  BP  tv  for  a  moderate  N 
is  often  substantially  improved  through  variance  reduction  techniques  such  as  importance  sampling, 
see,  e.g.,  [26].  If  all  realizations  v^, ...,  result  in  satisfactory  structural  performance  for  relevant 
designs,  i.e.,  g{x,  v-^)  <  0  for  all  j  =  1,  2, ...,  N,  then  globally  optimal  solutions  for  zq,  zi, ...,  zn  are 
all  zero.  This  implies  that  the  optimal  design  in  BP  tv  is  simply  the  x  that  minimizes  the  objective 
function  /(x)  over  X.  Consequently,  the  possibility  of  failure  of  the  structure  is  not  accounted  for 
in  BPtv  in  the  case  of  such  realizations.  Hence,  it  is  important  that  some  of  the  realizations  result 
in  g{x,v^)  >  0  for  relevant  designs.  We  can  typically  accomplish  this  by  increasing  N  or,  more 
efficiently,  by  importance  sampling,  which  we  describe  next. 

Let  W  be  a  random  vector  with  m  random  variables  with  joint  probability  density  function 
/w(w)  with  /w(w)  >  0  for  all  w  satisfying  /v(w)  >  0.  Let  w^,...,w'^  be  realizations  of  W. 
Then,  we  redefine 


BP/v  : 


min  /(x) 

X,Z 


S.t.  Zo  -I- 


N 

E 


(5r(x,  W^)  -  Zo) 


/v(w^) 

/w(wJ) 


<  0 

<  Zj,  j  =  l,2,...,N 
>  0,  j  =  l,2,...,N 


X  G  V. 


(23) 


By  generating  realizations  from  an  appropriately  selected  probability  density  /w(’w),  we  can  ensure 
that  a  substantial  number  of  realizations  w-^  satisfies  g{x,'w^)  >  0  for  relevant  designs.  In  practice, 
/w(w)  can  typically  be  selected  by  increasing  (decreasing)  mean  values  of  random  variables  de¬ 
scribing  loads  (material  strength).  For  more  sophisticated  approaches  to  selecting  /w(w)  we  refer 
to  [26]. 
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4  System  Reliability 

4.1  Problem  Formulation 


The  performance  of  a  structure  is  often  given  by  multiple  limit-state  functions  representing  quan¬ 
tities  such  as  stresses  and  deformations  at  different  locations.  Let  v),  k  =  1,2,...,  K,  be  a 
collection  of  limit-state  functions  describing  the  relevant  limit  states  for  a  structure.  We  define 
a  cut-set  to  be  a  (sub)set  of  these  limit-state  functions  with  the  characteristics  that  if  all  the 
limit-state  functions  in  the  cut-set  are  unsatisfactory  for  a  given  design  x  and  realization  v,  i.e., 
(/fc(x,  v)  >  0,  then  the  structure  experiences  system  failure.  A  cut-set  is  minimal  if  no  limit-state 
function  can  be  removed  from  the  cut-set  without  rendering  the  resulting  set  not  a  cut-set.  We 
refer  to  an  individual  limit-state  function  being  unsatisfactory  as  component  failure.  Suppose  there 
are  ic  minimal  cut-sets.  We  denote  the  set  of  limit-state  functions  belonging  to  minimal  cut-set  i 
by  Ci,  i  G  I  =  {1,2, ...,  ic}.  As  system  failure  occurs  in  the  event  of  component  failure  with  respect 
to  all  limit-state  functions  in  any  minimal  cut-set,  the  system  failure  probability  is  defined  as 


/(x)  =  P 


u  nte(x,v)>o} 


(24) 


If  the  cardinality  of  Ci,  denoted  \Ci\,  is  one  for  all  i  G  I,  then  the  structure  is  a  series  structural 
system  as  the  failure  of  any  component  results  in  system  failure.  On  the  other  hand,  if  ic  =  1,  then 
the  structure  is  a  parallel  system  as  system  failure  only  occurs  if  all  components  fail. 

It  follows  directly  from  (24)  that 


p®(x)  =  P[{5((x,  V)  >  0}], 


(25) 


where 


g{x,  v)  =  max  min  gk{x,  v 

k^Ci 


(26) 


is  a  system  limit-state  function.  Hence,  the  design  optimization  problem  with  system  failure  con¬ 
straints  generalizes  P  and  takes  the  form 


P®  :  min/(x) 

s.t.  p*(x)  <  1  —  ao 

X  G  A. 
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As  P®  is  at  least  as  intractable  as  P,  we  consider  a  formulation  involving  the  buffered  system  failure 
probability. 

4.2  Using  the  Buffered  Failure  Probability 

Following  the  approach  of  Section  3,  we  define  analogously  to  BP jv  the  problem 


where  we  use  in  (27)  the  fact  that 

maxmin5(fc(x,  v^)  -  zo  <  2:^,  j  =  l,2,  (28) 

iGl  keCi 

is  equivalent  to  (27).  The  relationship  between  P®  and  BP^  is  identical  to  those  between  P  and 
BPtv.  Hence,  we  recommend  designers  to  consider  BP^  instead  of  P®. 

In  the  case  of  series  structural  systems,  i.e.,  each  minimal  cut-set  includes  only  one  limit-state 
function,  BP^  is  identical  to  BP^r,  except  it  includes  more  constraints  of  the  same  form.  Hence, 
BP^  is  tractable  by  standard  nonlinear  optimization  algorithms  when  gk{^,  v),  k  =  1,2, ...,  K,  are 
continuously  differentiable.  Moreover,  convexity  is  preserved  as  in  the  case  of  BP^r. 

Cases  with  general  or  parallel  structural  systems  are  more  complicated.  The  minimum  over 
limit-state  functions  in  (27)  causes  BP^  to  become  a  nonsmooth  optimization  problem  even  if 
gk{x,v),  k  =  1,2,  ...,K,  are  continuously  differentiable.  Hence,  standard  nonlinear  optimization 
algorithms  are  not  applicable.  We  propose  three  alternative  approaches  to  overcome  this  difficulty. 

The  first  alternative  transcribes  the  problem  into  a  finite,  but  potentially  large  number  of 
optimization  subproblems.  Specifically,  BP^  is  equivalent  to 

min  min  f(x) 


28 


N 


S.t.  Zq  + 


^(1  -  ao) 

5fcy(x,v^)-zo  <  Zj,  i  e  I,j  =  1,2,...,N 
Zj  >  0,  j  =  l,2,...,N 
X  e  X. 


(29) 


This  problem  amounts  to  minimizing  Ilie/ 1^*1^  subproblems  essentially  of  the  form  BP^r  and 
retaining  the  design  with  the  best  objective  function  value.  A  main  advantage  of  this  transcription 
is  that  it  preserves  convexity.  That  is,  if  ^^(x,  v),  A:  =  1,  2,  are  convex  functions  with  respect 

to  X,  /(x)  is  a  convex  function,  and  X  is  a  convex  set,  then  each  of  the  Hie/  \Ci\^  subproblems 
are  convex. 

A  design  found  in  one  of  the  subproblems  can  be  used  to  warm  start  the  calculations  of  the  next 
subproblem.  However,  the  main  challenge  with  this  approach  is  the  large  number  of  subproblems 
to  solve.  If  it  is  not  practical  to  (approximately)  solve  all  subproblems,  then  it  is  always  possible  to 
solve  only  a  subset  of  the  subproblems.  This  provides  a  conservative  design  as  (29)  is  a  restriction 
of  (27)  and  further  improvement  might  be  possible  after  solving  other  subproblems. 

The  second  alternative  avoids  the  large  number  of  subproblems  by  using  exponential  smoothing 
[1,  17].  This  alternative  replaces  the  nonsmooth  function 

V^(x,  v)  =  min5rfc(x,  v).  (30) 

keCi 

in  (27)  by  a  continuously  differentiable  approximation.  For  any  approximation  parameter  e  >  0, 
let 

ffi(x,v;e)  =  -eln  I  ^  j  (31) 

\keCi  J 

be  this  approximation.  We  know  that 


0  < '!/'(x,v) -5i(x,v;e)  <  elnlCil  (32) 

for  all  X,  V,  and  e  >  0.  Hence,  the  smooth  approximation  5j(x,  V'^;e)  underestimates  '0(x,  v-^)  and 
the  error  in  the  approximation  vanishes  as  e  ^  0. 
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We  now  simply  replace  min^gc'^  gk{x,v^)  in  (27)  by  its  smooth  approximation  for  all  i  and  j. 
This  results  in  the  following  problem 


BP^(e)  ;  min/(x) 


s.t.  zo  + 


1 


N 


5j(x,  v^;e)  +  elnlCil  -  Zo  < 


ZJ  > 


X  G 


0 

Zj,  i  e  I,j  = 

0,  j  =  l,2,...,iV 

X. 


(33) 


Since  we  included  the  error  term  eln  IC*!  in  (33),  a  design  that  is  feasible  in  BP^(e)  is  also  feasible 
in  BP^.  If  the  limit-state  functions  (^^(x,  v),  k  =  1,2,...,  K,  are  continuously  differentiable,  then 
standard  nonlinear  optimization  algorithms  are  applicable  to  BP|^(e). 

We  observe,  however,  that  even  if  gki^,  v),  /c  =  1,  2,  ...,K,  are  convex,  BPf^(e)  is  not  a  convex 
optimization  problem.  In  essence,  minimal  cut-sets  with  cardinality  larger  than  one  introduce 
nonconvexity  in  the  design  optimization  problem.  We  also  note  that  exponential  smoothing  can  be 
used  in  BPjv  to  replace  the  N  constraints  (21)  by  one  single  constraints. 

The  third  alternative  for  solving  BP^  adapts  the  approach  in  [12].  In  that  paper  it  is  shown 
that  BP^  is  equivalent  to  the  following  problem: 

EBP^  :  min  /(x) 

x,z,/ih  ,i£l  ,j=l,...,N,k£Ci 


S.t. 


ZQ  + 


1 

iV(l  -  ao) 


N 


k'ij9kix,v^)  -  Zo 

keCi 


<  0 

<  Zj,  i  e  I,j  =  l,2,...,N 
>  0,  j  =  l,2,...,iV 


X  G  X 


(34) 


^  =  1,  i  £  I,j  =  1,...,N 

keCi 

k'ij  >  0,  i  G  I,  j  =  1, ...,  N,  k  G  Ci, 


where  gij,i  G  I,  j  =  1,  ...,iV,  k  G  Ci,  is  a  set  of  auxiliary  design  variables  that  effectively  “select” 
which  limit-state  functions  in  (27)  are  active.  The  equivalence  between  EBP|^  and  BP|^  is  in 
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the  sense  that  a  globally  (locally)  optimal  solution  from  one  problem  can  be  used  to  construct  a 
globally  (locally)  optimal  solution  of  the  other  problem.  If  gki^,  v),  /c  =  1,  2, K,  are  continuously 
differentiable,  then  standard  nonlinear  optimization  algorithms  are  applicable  for  solving  EBP^. 
However,  even  if  gk{x,  v),  A:  =  1,  2, ...,  K,  are  convex,  EBP^  may  not  be  a  convex  problem  because 
(34)  involves  a  product  of  design  variables. 

5  Computational  Studies 

We  illustrate  the  use  of  the  buffered  failure  probability  with  the  design  of  a  truss  structure  and  a 
motor  vehicle. 

5.1  Optimal  Truss  Design 

Consider  the  simply  supported  truss  in  Figure  11.  Let  14  be  the  yield  stress  of  member  k,  k  = 
1,  2, ...,  7.  Members  1  and  2  have  lognormally  distributed  yield  stresses  with  mean  100  N/mm^  and 
standard  deviation  20  N/mm^.  The  other  members  have  lognormally  distributed  yield  stresses  with 
mean  200  N/mm^  and  standard  deviation  40  N/mm^.  The  yield  stresses  of  members  1  and  2  are 
correlated  with  correlation  coefficients  0.8.  However,  their  correlation  coefficients  with  the  other 
yield  stresses  are  0.5.  Similarly,  the  yield  stresses  of  members  3-7  are  correlated  with  correlation 
coefficients  0.8.  The  truss  is  subject  to  a  random  load  Vs  in  its  mid-span.  Vs  is  lognormally 
distributed  with  mean  1000  kN  and  standard  deviation  400  kN.  The  load  Vg  is  independent  of  the 
yield  stresses.  We  use  a  joint  lognormal  distribution  (see  [5],  Section  7.2)  and  the  above  correlation 
coefficients  to  approximate  the  joint  distribution  of  V  =  (Vi,  14,  •••,  Lg). 

The  design  vector  x  =  (xi,  X2, ...,  xy),  where  Xk  is  the  cross-section  area  (in  1000  mm^)  of 
member  k.  The  truss  fails  if  any  of  the  members  exceed  their  yield  stress.  (We  ignore  the  possibility 
of  buckling.)  This  gives  rise  to  seven  limit  state  functions; 

5'fc(x,v)  =  vs/Ck-VkXk,k  =  1,2,. ..,7,  (35) 

where  Cfc  is  a  factor  given  by  the  geometry  and  loading  of  the  truss.  From  Figure  11,  we  determine 
that  Ck  =  l/(2\/3)  for  k  =  1,2,  and  (k  =  l/\/3  for  k  =  3,4, ...,  7. 
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Figure  11:  Design  of  Truss 


We  impose  the  constraint  that  the  series  system  failure  probability  with  the  seven  limit-state 
functions  should  be  no  larger  than  0.00135,  i.e., 


p^(x)  =  P 


U{fffc(x,V)  >0} 


<  0.00135. 


(36) 


LA:=1  J 

We  also  impose  the  14  deterministic  constraints  0.5  <  <  2,k  =  1,2,...,  7,  that  limit  the  allowable 

area  of  each  member  to  be  between  500  mm^  and  2000  mm^.  We  seek  a  design  of  the  truss  that 
minimizes  the  cost  of  the  truss.  Since  all  members  are  equally  long,  the  cost  is  /(x)  =  J2k=i^k- 
This  problem  is  of  the  form  and,  hence,  we  solve  BP^  with  (27)  replaced  by 


9fc(x,  V-')  -zo<  Zj,  k  =  l,2,  ...,7,j  =  1,2,  ...,iV 


(37) 


as  our  example  is  a  series  structural  system.  Since  the  limit-state  functions,  objective  functions, 
and  constraints  are  linear  in  x,  BP^  is  a  linear  program  that  can  be  solved  quickly  by  standard 
optimization  solvers.  We  use  sample  size  N  =  10000  and  the  variance  reduction  technique  of 
Subsection  3.4  where  we  select  the  sampling  distribution  to  be  identical  to  the  original  distribution 
except  that  the  mean  value  of  the  load  is  increased  with  three  standard  deviations. 


32 


We  implement  the  resulting  linear  instance  of  BP^  in  the  General  Algebraic  Modeling  System 
(GAMS)  Distribution  22.9  [6]  on  a  laptop  computer  with  1  GB  of  RAM  and  2.16  GHz  processor 
running  Windows  XP.  The  globally  optimal  design  of  is  found  by  the  solver  CPLEX  11.2  [11] 
with  default  options  in  19.4  seconds.  The  optimized  design  is  shown  in  row  3  of  Table  1  with  the 
resulting  buffered  failure  probability  in  the  last  column  of  the  table.  In  this  case,  the  buffered  failure 
probability  is  about  three  times  larger  than  the  failure  probability.  Here,  the  failure  probability  is 
estimated  by  importance  sampling  with  a  5%  coefficient  of  variation  using  an  independent  sample. 
For  a  rigorous  solution  validation  procedure  we  refer  to  [21].  For  comparison,  we  also  report  in 
row  4  of  Table  1  the  design  found  for  the  same  truss  in  [25]  by  approximately  solving  P®  using 
sample  average  approximations.  Since  P®  requires  the  failure  probability  to  be  no  larger  than  a 
threshold  and  BP^,  effectively,  imposes  the  same  threshold  on  the  buffered  failure  probability,  the 
design  of  row  4  is  naturally  cheaper  than  the  one  in  row  3.  However,  the  former  design  is  less  safe 
with  an  estimated  failure  probability  of  0.00153  (5%  coefficient  of  variation  of  estimate),  which 
slightly  exceeds  the  threshold  of  0.00135.  While  the  algorithm  in  [25]  is  guaranteed  to  converge 
to  a  feasible  design  satisfying  the  KKT  conditions  under  suitable  assumptions,  termination  of  the 
algorithm  after  a  finite  amount  of  calculation  time  may  result  in  such  infeasibilities.  In  contrast, 
the  design  obtained  by  using  the  buffered  failure  probability  has  an  estimated  failure  probability 
below  the  required  threshold  of  0.00135.  Moreover,  the  calculation  time  of  the  algorithm  in  [25]  is 
substantially  longer  than  that  of  solving  BP|^,  with  a  time  exceeding  one  hour  to  obtain  the  design 
in  row  4  of  Table  1.  While  an  improved  implementation  of  the  algorithm  in  [25]  will  reduce  this 
time,  the  advantage  of  the  buffered  failure  probability  appears  substantial. 

To  better  compare  the  design  obtained  using  the  buffered  failure  probability  with  that  using  the 
failure  probability,  we  also  solve  BP^  with  probability  threshold  0.00410.  This  threshold  equals 
the  buffered  failure  probability  of  the  design  in  row  4  of  Table  1.  Row  5  of  the  table  gives  the 
resulting  design  obtained  after  20.5  seconds  using  CPLEX.  We  see  that  the  designs  in  rows  4  and  5 
are  essentially  identical,  which  indicate  that  optimization  with  the  buffered  failure  probability  gives 
a  similar  design  to  that  obtained  using  the  failure  probability  when  the  threshold  is  appropriately 
adjusted.  We  note  again  that  the  computing  time  is  dramatically  reduces  when  using  the  buffered 
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Method 

1  —  ao 

1 

Design  of  member 
2  3  4 

(in  mm^) 

5  6 

7 

Cost 

(mm^) 

Failure 

Prob. 

Buffered 

Prob. 

BP^ 

0.00135 

1320 

1332 

1272 

1278 

1271 

1278 

1271 

9022 

0.00047 

0.00135 

[25] 

0.00135 

1138 

1156 

1118 

1107 

1119 

1113 

1108 

7859 

0.00153 

0.00410 

BP^ 

0.00410 

1153 

1179 

1100 

1105 

1106 

1109 

1101 

7852 

0.00154 

0.00410 

Table  1:  Design  of  truss 


failure  probability. 

5.2  Motor  Vehicle  Design 

We  consider  an  example  given  in  [27]  (see  also  [8])  where  the  goal  is  to  minimize  the  weight  of  a 
part  of  a  motor  vehicle  subject  to  reliability  constraints  related  to  side  impact.  We  formulate  this 
problem  in  the  form  P®  with  a  series  system  failure  probability  with  respect  to  ten  limit-states 
functions  and  a  reliability  level  1  —  ao  of  0.0013.  The  limit-state  functions  are  surrogate  models 
of  the  real  structural  performance;  see  [27].  The  example  has  seven  design  variables  relating  to 
the  thickness  of  material.  (The  paper  [27]  includes  four  additional  variables,  which  we  simply  fix 
to  the  values  reported  in  [27],  i.e.,  0.345,  0.345,  0,  and  0.)  All  thicknesses  must  be  in  the  interval 
[0.5  1.5].  The  thicknesses  cannot  be  manufactured  exactly  and,  hence,  the  limit-state  functions 
include  normally  distributed  manufacturing  errors  with  zero  mean  and  standard  deviation  0.03  for 
each  thickness.  The  errors  are  statistically  independent.  We  refer  to  [27]  for  details  of  this  example. 

We  implement  BP|^  for  this  example  with  sample  size  N  =  7500  using  the  same  hardware 
as  above,  but  now  solve  the  problem  using  SNOPT  [9]  as  implemented  in  TOMLAB  [30].  Table 
2  gives  the  optimized  design  in  row  3,  which  was  obtained  after  166  seconds,  and  the  resulting 
buffered  failure  probability;  see  the  last  column.  The  corresponding  failure  probability  is  estimated 
by  Monte  Carlo  sampling  with  a  5%  coefficient  of  variation  using  an  independent  sample;  see  the 
second  to  last  column  of  row  3.  For  comparison,  we  also  report  the  design  given  in  [27]  with 
estimated  failure  and  buffered  failure  probabilities  (5%  coefficient  of  variation);  see  row  4  of  Table 
2.  Again,  we  see  that  our  methodology  results  in  a  reasonable  design  in  short  computing  time. 
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Method 

Optimized  Design 

Cost 

Failure 

Buffered 

Xl 

X2 

X3 

X4 

X5 

X6 

X7 

Prob. 

Prob. 

BP^ 

0.5000 

1.3524 

0.5000 

1.2989 

0.6103 

1.5000 

0.5000 

24.60 

0.00067 

0.00130 

[27] 

0.5000 

1.3251 

0.5000 

1.2919 

0.5964 

1.5000 

0.5000 

24.37 

0.00347 

0.01769 

Table  2;  Design  of  Motor  Vehicle 


6  Conclusions 

We  discuss  several  theoretical,  practical,  and  computational  issues  associated  with  the  failure  prob¬ 
ability  with  particular  emphasis  on  difficulties  arising  in  design  optimization.  We  propose  an  al¬ 
ternative  measure,  the  buffered  failure  probability,  that  offers  significant  advantages.  The  buffered 
failure  probability  accounts  for  the  degree  of  violation  of  a  performance  threshold,  is  more  conserva¬ 
tive  than  the  failure  probability,  and  is  handled  with  relative  ease  in  design  optimization  problems. 
The  paper  formulates  several  design  optimization  problems  in  terms  of  the  buffered  failure  probabil¬ 
ity  and  discusses  their  relation  to  design  optimization  problems  in  terms  of  the  failure  probability. 
We  find  the  buffered  failure  probability  to  be  superior  to  the  failure  probability  and  recommends 
its  use  in  design  and  optimization  of  structures. 

While  the  buffered  failure  probability  appears  promising  for  use  in  design  optimization  with 
reliability  constraints,  its  applicability  in  other  optimization  models  such  as  those  with  a  von 
Neumann-Morgenstern  maximum  expected  utility  criterion  is  unclear.  Moreover,  the  buffered  fail¬ 
ure  probability  requires  the  estimation  of  an  expectation,  which  may  be  computationally  costly,  and 
may  result  in  large-scale  optimization  models.  These  challenges  should  be  the  subject  of  further 
study. 
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